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We use extensive first-principles quantum mechanical calculations to show that, although the 
static lattice and harmonic vibrational energies are almost identical, the anharmonic vibrational 
energy of hexagonal ice is significantly lower than that of cubic ice. This difference in anharmonicity 
is crucial, stabilising hexagonal ice compared with cubic ice by at least 1.4 meV/H20, in agreement 
with experimental estimates. The difference in anharmonicity arises predominantly from molecular 
0-H bond stretching vibrational modes and is related to the different stacking of atomic layers. 


I. INTRODUCTION 

Six-fold symmetric snow crystals are formed from 
hexagonal ice (Ih), which covers about 10% of the Earth’s 
surface and plays a prominent role in determining its cli¬ 
mate M. A cubic form of ice also occurs in nature, 
but is very rare ii- The structures of pure hexagonal 
and cubic ice differ only in the stacking of layers of tetra- 
hedrally coordinated water molecules (see Fig. [1]). Yet, 
hexagonal ice is thermodynamically more stable than cu¬ 
bic ice, with experiments indicating the difference in sta¬ 
bility to lie in the meV/H 20 range There is a 

growing realisation that real “cubic ice” typically con¬ 
tains many stacking faults and is not pure cubic ice (Ic) 
as originally suggested by Konig [TgI . Stacking faulted ice 
is a highly complex material, whose nature and proper¬ 
ties depend heavily on the free energy difference between 
pure Ih and Ic. 

The very similar free energies of Ih and Ic have so 
far prevented state-of-the-art first-principles quantum 
mechanical calculations from explaining the stability of 
Ih. Both density functional theory (DFT) and diffusion 
quantum Monte Carlo studies have found that Ih and Ic 
are almost degenerate in energy when nuclear motion is 
neglected El- Our calculations show that the harmonic 
zero-point vibrational energies of Ih and Ic are large, 
at roughly 700meV/H20, but they are almost identi¬ 
cal. Consequently, when averaged over different proton- 
orderings, the two phases are almost degenerate when 
harmonic vibrations are included (see Fig. ED- However, 
the small mass of hydrogen gives rise to large amplitude 
vibrations and large anharmonic effects. 

A substantial body of theoretical work exists on wa¬ 
ter and ice, based on force-field path-integral molecular 
dynamics (PIMD), first-principles classical molecular dy¬ 
namics (MD), and first-principles vibrational calculations 
in the quasi-harmonic approximation. This work has led 
to significant successes in understanding the important 
role of quantum vibrations and anharmonicity for vari¬ 
ous phenomena observed experimentally. Examples in¬ 
clude, but are not limited to, (a) the isotope effects, e.g., 
on the melting temperature, of Ih upon going from pro- 
tiated to deuterated ice [iilii, (b) accurate 0-H bond 
lengths and infrared 0-H stretching frequencies in water 


(c) reproduction of the anomalous thermal ex¬ 
pansion, and the isotope effect on the volume in Ih 
and (d) the heat capacity of water [25| . Attempts to cal¬ 
culate the relative stability of Ih and Ic have either relied 
on empirical force-fields such as TIP4P [26|, or have 
lacked an accurate description of anharmonicity mil. 
TIP4P has since been shown to produce incorrect proton¬ 
ordering energetics and an incorrect static lattice energy 
difference between Ih and Ic com par ed to highly accurate 
diffusion Monte Carlo methods Moreover, no suc¬ 

cessful attempts to explain the origin of the greater sta¬ 
bility of Ih have been made. With our fully anharmonic, 
first-principles DFT study we show that the inclusion of 
accurate anharmonic quantum nuclear motion is decisive 
in stabilising Ih with respect to Ic, and relate the differ¬ 
ence in stability to the different stacking of the atomic 
layers. 

The water molecules in both Ih and Ic are tetrahedrally 
coordinated, each donating and accepting two hydrogen 
bonds and thus satisfying the “Bernal-Fowler ice rules” 
[Slj . The oxygen sublattices of Ih and Ic arise from an 
ABAB and an ABC stacking of puckered layers of oxygen 
atoms, respectively. Correspondingly, the basic building 
blocks of Ih are chair- and boat-form hexamers, while 
Ic is built exclusively from chair-form hexamers (Fig. 
[T]). Pure Ic has so far proven elusive. Experimentally 
synthesised “Ic” typically contains many stacking faults 
which strongly affect its physical and chemical proper¬ 
ties i m [ 3 ^ . Ice containing cubic sequences inter¬ 
laced with hexagonal sequences is commonly known as 
stacking-disordered ice (Isd). Unlike Ih and Ic, which re¬ 
fer to a unique stacking arrangement of puckered layers, 
Isd refers to the infinite set of possible stacking sequences. 
This set smoothly connects Ih as one end member to Ic as 
the other. Isd has trigonal P3ml symmetry jssMsbj. As 
of yet, it is unclear whether stacking disorder is kineti- 
cally or thermodynamically driven. A full understanding 
of Isd will require understanding the properties of Ih and 
Ic, including their relative stability. At the most basic 
level, the free energy difference between Ih and Ic is re¬ 
quired to understand why Isd is found to anneal to Ih 
rather than Ic. 

The fraction of cubic stackings of layers, or “cubic- 
ity”, typically does not exceed around 60% in exper¬ 
iments and both the fraction itself as well as the nature 
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FIG. 1. The oxygen sublattices in (a) hexagonal ice, Ih, and 
(b) cubic ice, Ic, consist of ABAB and ABC stacked puckered 
layers (red, blue and green colours), respectively. Projections 
of the periodic structures along the direction orthogonal to the 
puckered layers are shown on the lower faces of the diagrams. 
Ih is built of the chair- and boat-form hexamers shown in a-I 
and a-II, respectively, while Ic is built exclusively from chair- 
form hexamers as in b-I. The structures were visualised using 
VESTA 0. 


of the stacking arrangements depend heavily on the syn¬ 
thesis pathway [Ij, . Crucially, ice synthesised via 

both homogeneous and heterogeneous freezing of (super¬ 
cooled) water has a random stacking of cubic and hexag¬ 
onal layers lallai, which is consistent with a layer-by- 
layer growth mechanism. Heterogeneous freezing in par¬ 
ticular is central to atmospheric and climate physics and, 
due to random stacking-disorder, depends vitally on the 
free energy difference between Ih and Ic. It is clear that 
Isd is an extremely important and highly complex mate¬ 
rial. 

Molecular dynamics and Monte Carlo simulations us¬ 
ing empirical ice potentials to model ice nucleation pro¬ 
cesses have successfully reproduced stacking-disorder (see 
[23,[36|,[33 and references therein). However, they strug¬ 
gle, amongst other issues, with the accuracy of the em¬ 
pirical potentials in describing the melting temperature 
and the relative stability of Ih and Ic [38[ . 

In the following we limit ourselves to the study of pure 
Ih and Ic. 


II. COMPUTATIONAL MODEL 

Atomistic simulations of proton-disordered systems 
such as Ih and Ic require sets of explicit atomic positions 
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FIG. 2. Harmonic (Char, empty squares, lower part) and an- 
harmonic (Ganh, filled squares, upper part) energies including 
zero-point nuclear motion for different Ih (blue) and Ic (red) 
proton-orderings (Supplementary Table 1). The averages over 
proton-orderings, Ghar and Ganh, are shown as thick horizon¬ 
tal dotted and solid lines, respectively. All energies are mea¬ 
sured with respect to G^iJ^ (Ih Gmc2i, structure 15). 


and a calculation must be performed for each proton¬ 
ordering studied. The number of proton-ordered, ener¬ 
getically quasi-degenerate structures allowed by the ice 
rules increases exponentially with the size of the simula¬ 
tion cell. This leads to Pauling’s residual configurational 
entropy 1391 . ^onfig, which has been confirmed experi¬ 
mentally [40, . For large systems with negligible sur¬ 
face effects the associated configurational free energies 
of Ih and Ic, AGconfig = —FS'config, are almost identical 
[ 421 , 143 . We therefore neglect AGconfig in the following. 

To gain an understanding of the effects of proton¬ 
ordering on the vibrational properties of ice, we consider 
16 distinct proton-ordered eight-molecule Ih configura¬ 
tions as constructed by Hirsch and Ojamae m , and 11 
distinct proton-ordered eight-molecule Ic configurations 
El- We also consider the “conventional” hexagonal, 12- 
molecule PGscm Ih and quasi-cubic, eight-molecule Pds 
Ic structures (numbers 13 and 1 in Fig. [21 respectively). 
Details of the proton-ordered structures and their num¬ 
bering are provided in Supplementary Table 1. 

We performed electronic structure calculations using 
plane-wave pseudopotential DFT as implemented in the 
CASTEP code [13 (version 7.02). We employed the 
Perdew-Burke-Ernzerhof (PBE) [^ generalised gradient 
approximation to the exchange-correlation functional, 
and on-the-fiy generated ultrasoft pseudopotentials [13 
with core radii of 0.7 A and 0.8 A for the hydrogen and 
oxygen atoms, respectively. Supplementary Section V 
describes results obtained with other density function¬ 
als. We used a plane-wave energy cut-off of 1600 eV 
and Monkhorst-Pack reciprocal space grids of spacing 
less than 27 r x 0.04 A for all total energy calculations 
and geometry optimisations. The resulting energy dif¬ 
ferences between frozen-phonon configurations are con¬ 
verged to below 10 “^eV/H 2 O, the atomic positions are 
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converged to within 10“^A, and the residual forces to 
within 10“^eV/A. 

We obtained harmonic vibrational free energies from 
the harmonic frequencies calculated using the k-space 
Fourier interpolated dynamical matrix. The latter was 
obtained by Fourier transforming the real-space matrix 
of force constants constructed using a finite displacement 
method. Anharmonic vibrational free energies were cal¬ 
culated using the method described in [^, which has 
so far been successfully applied to high-pressure systems 
[H, [ 5 ^. As in [m, we describe the SAT-dimensional BO 
energy surface (where M is the number of atoms in the 
simulation cell) by mapping ID subspaces along the har¬ 
monic normal mode axes up to large amplitudes of four 
times the harmonic root-mean-square (RMS) displace¬ 
ments, where anharmonicity is important. We then re¬ 
construct the SA/"-dimensional BO surface from the ID 
subspaces. The resultant representation of the BO energy 
surface is an approximation to the true SA/"-dimensional 
BO energy surface. This approximation only weakly af¬ 
fects the free energy difference between Ih and Ic (see 
Supplementary Section VI). The ID energy surfaces were 
fitted using cubic splines. The anharmonic vibrational 
Schrodinger equation was solved within a vibrational self- 
consistent field (VSCF) framework. The vibrational wave 
function was expanded in a basis of simple harmonic os¬ 
cillator eigenstates, and the inclusion of 25 states for each 
vibrational degree of freedom was found sufficient to ob¬ 
tain converged results. 


III. RESULTS 


Our calculations show that the static lattice energies of 
Ih and Ic, F’static, vary by up to 5 meV/H 2 O with proton¬ 
ordering. This agrees with Refs. [13, H, and, more 
importantly, with Ref. [52|, which evaluates DFT static 
lattice energies for 16 8-molecule orthorhombic, 14 12- 
molecule hexagonal and 63 48-molecule orthorhombic Ih 
proton-orderings. This strongly suggests that our sets 
of proton-orderings provide a good representation of the 
distribution of energies in disordered ice. 

We also find that the harmonic vibrational contribu¬ 
tions to the free energies of different proton-orderings, 
AGhar, vary by up to 2 meV/H 2 O. We have employed the 
VSCF method described above to calculate the anhar¬ 
monic contribution to the vibrational free energy, AGanh, 
finding a variation between proton-orderings of up to 
5 meV/H2 0 

The free energies of Ih and Ic at the harmonic vibra¬ 
tional level, Ghar = ^static + AGhar, averaged over the 
proton-orderings, are virtually indistinguishable: 


^ - GI*;, = 0.2 ± 2.4meV/H20, 

where the quoted errors are the RMS variations across 
different proton-orderings. The anharmonic energies of 
the Ih configurations, AGanh, on the other hand, though 
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FIG. 3. Difference in the anharmonicity of the BO energy 
surface, AFanh, of Ih Cmc2i and Ic lAimd for various func¬ 
tionals (thin solid and dotted lines). The dominant contribu¬ 
tion to the anharmonicity is quart ic. The difference between 
the results from different functionals is so small it can barely 
be distinguished. The different widths of the vibrational den¬ 
sities, |Tanh|^, of Ih and Ic (blue and red thick dotted lines, 
respectively) is indistinguishable on the scale of the figure. 



also positive, are systematically lower than those of the 
Ic configurations, so that the total free energy, Ganh = 
^static + AGhar + AGanh, averaged over the different 
proton-orderings, is significantly lower for Ih than for Ic: 

Air'^Ganh ^ ^ ^ = 6.5 ± 3.1 meV/H20. 

The values obtained for A^^“^^^Ghar and A^^“^^^Ganh de¬ 
pend significantly on the method of averaging. For exam¬ 
ple, using a Boltzmann distribution for the free energy of 
the proton-orderings leads to values of Ag^^mann^har ~ 
5.5meV/H20 and 6.1 meV/H20 at 10 and 100 K, respec¬ 
tively. 

It is noteworthy that, given cell volumes that are rea¬ 
sonably close to experiment, the differences in Fgtatic, 
AGhar and AGanh between Ih and Ic depend only weakly 
on the details of the DFT calculations and in particular 
on the choice of exchange-correlation functional (see Fig. 
[ 3 ] and Supplementary Section V). 

The most stable proton-ordered configurations of Ih 
and Ic, referred to as Xlh and XIc, display free energy 
differences very similar to A^^“^^^Ghar and A^^“^^^Ganh- 
However, the inclusion of anharmonic vibrational ener¬ 
gies changes the relative stability of the different proton- 
orderings. Xlh is experimentally and theoretically known 
to have space group Cmc2i [3, (structure 15 in Fig. 
ED, a result confirmed by our calculations. A structure 
of space group Mimd has been proposed 0 (struc¬ 
ture 2) for XIc on theoretical grounds. Our calculations 
support this proposal at the harmonic vibrational level, 
but the inclusion of anharmonic contributions suggests 
that Ic Pc, Pca2i and Pna2i (structures 7, 9 and 12) 
may also be strong candidates for XIc. Ref. [HJ reports 
experimental evidence for Ic lAimd and Pna2i in par¬ 
tially proton-ordered Ic via Fourier transform infrared 
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spectroscopy. This lends significant support to our re¬ 
sult that anharmonicity provides the decisive contribu¬ 
tion to the energy differences between proton-orderings, 
since Pna2i has a high free energy at the static lattice 
and harmonic vibrational levels and only becomes a low- 
free-energy structure when anharmonicity is taken into 
account. 

At typical experimental temperatures of below 100 K, 
the proton-ordering is largely frozen in. Consequently 
one cannot expect to measure a change in free energy cor¬ 
responding to a transition from Ic to Ih in thermal equi¬ 
librium, but rather (assuming the Ic sample is annealed 
at low temperatures and consists mostly of XIc) a change 
in free energy corresponding to transitions from XIc to a 
proton-ordering of Ih, which is likely to be smaller than 
As indicated in Fig. [21 we evaluate a lower 
bound on the free energy difference as 

Amht^‘’G'anh = min {G^'=„h} - max 
= 1.4±0.3meV/H2O. 

This lower bound is consistent with, but on the high 
side of, experimentally measured free energy differences 
of 0.3- l.GmeV/HsO Notably, the experimental 

value is rather uncertain, mainly because the free en¬ 
ergy difference is very small, and because the Ic samples 
are typically not fully characterised in terms of stacking 
faults or proton ordering. 



FIG. 4. Ganh for H 2 O as a function of temperature. The solid 
blue and red lines indicate results for Ih PGscm and Ic P43 
simulation cells containing 96 and 64 molecules, respectively. 
The light blue and light red areas indicate the spread in free 
energies arising from the different Ih and Ic proton-orderings. 

As shown in Fig. (H this energy difference remains 
roughly constant from zero temperature up to 273 K and 
thus stabilises Ih over a wide range of temperatures. 

In D 2 O, the heavier mass of deuterium firstly leads 
to reduced vibrational frequencies of the harmonic vi¬ 
brational modes and thus reduced vibrational energies. 
For Xlh D 2 O, AGhar = 507.67 meV/D20 compared 
to 692.34 meV/H20 for Xlh H 2 O. Secondly, D 2 O has 
smaller vibrational amplitudes and consequently smaller 



FIG. 5. Zero temperature (empty squares, lower part) 

and (filled squares, upper part) for different deuterated 

Ih (blue) and Ic (red) proton-orderings. The averages over 

proton-orderings, G^^^ and are shown as thick dotted 

and solid lines, respectively. All energies are measured with 
respect to G^^^ of Ih Gmc2i. 


anharmonic vibrational energies than H 2 O, resulting in 
a smaller difference between the free energies of Ih and 
Ic, 


A Ic^Ih^D20 _ ^Ic-D 20 ^Ih-D20 

^av '-^anh — '-^anh '-^anh 


3.7±2.9meV/D20, 


as shown in Fig. O This suggests that it could be easier 
to synthesise deuterated than protiated cubic ice. More¬ 
over, the smaller anharmonic energy calculated for D 2 O 
compared with H 2 O implies that the most likely candi¬ 
date for the ground state of cubic D 2 O has lAimd sym¬ 
metry (structure 2) and is thus different from the ground 
state of H 2 O. The predicted ground state of hexagonal 
D 2 O is Gmc2i, as for H 2 O. 

We have investigated the convergence of AGhar and 
AGanh with respect to the size of the simulation cell using 
cells containing up to 192 molecules for Ih PGscm and 128 
molecules for Ic Pds as shown in Fig. |6l 

The vibrational energies of the different proton- 
orderings were calculated using 64-molecule cells for the 
harmonic vibrational energy and eight-molecule cells for 
the anharmonic energies. The latter is justified by the 
short-range nature of anharmonicity, evidence for which 
is described in section HYi Calculations using cells with 
up to 192 molecules for the PGscm Ih structure and 128 
molecules for the Pds Ic structure indicate that the differ¬ 
ence between AGhar for Ih and Ic is converged to within 
0.1meV/H2O using 64-molecule simulations cells. The 
anharmonic corrections converge analogously, and the 
difference between AGanh for the two structures is con¬ 
verged to within 0.2 meV/H 2 O using 64-molecule simula¬ 
tions cells. More details may be found in Supplementary 
Section II. 
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FIG. 6. The harmonic free energy, AGhar, and anharmonic 
free energy, AGanh, are shown as functions of the simulation 
cell size in the upper and lower panels, respectively. 


IV. DISCUSSION 


The origin of the difference in the anharmonicities of Ih 
and Ic can be traced back to the libration (80 — 140 meV) 
and, predominantly, the molecular symmetric and anti¬ 
symmetric 0-H bond stretching modes (365 — 380 meV 
and 395 —410 meV, respectively) indicated in Fig. [71 The 
vibrational density of states (DoS) and the distribution of 
anharmonic corrections over the vibrational frequencies 
show (Fig. IH]) that the dominant anharmonic contribu¬ 
tions arise from the 0-H bond stretching modes, which 
correspond to large amplitude displacements of the hy¬ 
drogen atoms relative to their neighbouring (essentially 
stationary) oxygen atoms. The comparatively small role 
of the oxygen atoms is confirmed by studying Ih and Ic 
analogues with fixed oxygen positions, which recover the 
value of AGanh observed for real Ih and Ic to within 5% 
(see Supplementary Table III). The 0-H bond stretch¬ 
ing modes contribute > 2/3 of the difference in anhar- 
monicity between Ih and Ic, when averaged over proton- 
orderings. Note that the energy difference between Ih 
and Ic shown in Fig. [8] increases at high frequencies. 

Variations in vibrational frequencies with cubicity have 
recently been identified in infrared absorption experi¬ 
ments [6[. Carr et al observed that the 0-H stretching 
modes were shifted to higher vibrational frequencies with 
increasing cubicity of their Isd samples. They also ob¬ 
served an increasing broadening of the absorption peak. 
According to Carr et a/., both trends are thought to be 
associated with the stacking disorder, which peaks at a 
cubicity of 50%. For samples with cubicities of 50% Carr 
et al. observed shifts of 28 ± 2 cm“^ and 13 ± 2 cm“^ for 
protiated and deuterated Isd, respectively. Our calcula¬ 
tions reproduce the widening of the 0-H stretching peak 


(a) Ih vibrational bandstructure and DoS. 
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(b) Ic vibrational bandstructure and DoS. 
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FIG. 7. Vibrational bandstructure and DoS of Ih (a) and Ic 
(b) split into crystal modes (I), libration modes (II), molecular 
0-H bending modes (III) and molecular 0-H anti-symmetric 
and symmetric stretching modes (IV and V, respectively). 
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in the vibrational DoS and produce shifts of 70 ± 5 cm“i^ 
and 30 ±5 cm“i^ for protiated and deuterated Isd, respec¬ 
tively. I.e., at about double the cubicity of Carr et a/., 
we calculate shifts that are about twice as large as those 
of Carr et al. Thus our values are consistent with the re¬ 
sults of Carr et al. if the shift scales roughly linearly with 
the degree of cubicity rather than the amount of stacking 
disorder, and is largest in pure Ic. In accurately repro¬ 
ducing the ratio between the blueshifts of the molecular 
stretching frequencies for the protiated and deuterated 
phases, our results also provide a good account of the 
effects of isotopic substitution. 

The role of the high energy modes can be further illu¬ 
minated by considering the H-H radial distribution func¬ 
tions (Fig. [9] (b)) of Ih and Ic and RMS displacements 
of the protons (see Supplementary Table IV). The lat¬ 
ter show that the harmonic vibrational amplitudes in Ic 
are about 1% smaller than in Ih. Yet, at the harmonic 
level the H-H radial distribution functions of Ih and Ic 
are essentially identical for distances less than about 3 A. 
Beyond 3 A the static structures of the two phases differ. 
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FIG. 8. Cumulative anharmonic energies and vibrational DoS 
of H 2 O averaged over proton-orderings as a function of fre¬ 
quency. The upper panel shows the anharmonic energies for 
Ih Cmc 2 i and Ic PAs as solid blue and red lines and the 
spread in free energies in Ih and Ic due to different proton- 
orderings as light red and light blue areas. The lower panel 
shows that the harmonic DoS of Ih and Ic (dotted lines) are 
practically identical. The larger increase in vibrational energy 
due to anharmonicity in Ic than in Ih is reffected in the shifts 
of their anharmonic DoS (continuous lines) with respect to 
the harmonic DoS. 
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FIG. 9. The RDFs for Ih PQscm (blue) and Ic PA 3 (red) 
are sampled from the anharmonic nuclear wavefunctions. For 
radii < SA the difference between Ih and Ic at the harmonic 
level (green), A^har = 5'har~5'har? minimal. The difference 
at the anharmonic level (black), A^anh = is small, 

but non-negligible. The features in A^har and A^anh for radii 
beyond around 3 A originate predominantly from differences 
in the static structures of Ih and Ic. 


Unlike the RMS displacements of the protons, the RDFs 
measure a two-particle quantity that describes correlated 
motions of pairs of protons. At the harmonic level the 
protons in Ic move less with respect to their equilibrium 
positions than in Ih, but they move by just as much with 
respect to each other. We note that the 0-H and H-H 
RDFs for Ih shown in Figs. [9] (a) and (b) agree well with, 
e.g., the experimental RDFs in [HBj . 

While the protons in both phases feel the same local 
environment, differences occur starting with the fourth- 
nearest-neighbour protons (see Fig. [ini). Also, for sys¬ 
tems as small as 8 to 12 molecules AGanh is already 
about 3/4 of the converged value (Supplementary Fig. 
1). This system size limits the wavelength of the vibra¬ 
tional modes responsible for the difference in anharmonic 
energies to roughly the same distance as the separation of 
fourth-nearest-neighbour protons. Together these obser¬ 
vations indicate that the influence of more distant nuclei 


is small. Allowing for both chair- and boat-form hexam- 
ers, there are 12 distinct arrangements of fourth-nearest- 
neighbour pairs of protons. Out of these 12 arrangements 
only 8 are realised in the proton-orderings we have con¬ 
sidered. These are shown in Fig. [TTl Three of these 
arrangements (numbers VI - VIII) are associated with 
boat-form hexamers of H 2 O molecules, which only ex¬ 
ist in Ih, and 5 (numbers I - V) are associated with the 
chair-form hexamers of H 2 O molecules found in Ih and 
Ic. The RDFs in Fig. [TT] show that, on going from Ic 
to Ih, arrangements II and IV are depopulated. For ar¬ 
rangements II and IV the displacement of the first proton 
of the pair from its equilibrium position along its hydro¬ 
gen bridge bond leads to a large displacement relative to 
the second proton. Conversely, arrangements VI - VIII, 
for which the same displacement of the first proton leads 
to a far smaller displacement relative to the second, are 
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(a) Typical chair-form configuration found in Ih and Ic. 



is typically driven by a strong quartic contribution to 
the respective BO energy surfaces (see Supplementary 
Fig. 6). Moreover, at the anharmonic level the RDFs 
show that the protons in ice Ic move less with respect to 
each other than in Ih, instead of just moving less with 
respect to their equilibrium positions, as they do at the 
harmonic level. This difference in the relative motion of 
pairs of protons is the origin of the difference in AGanh 
between Ih and Ic. The larger effect of anharmonicity 
in Ic is again due to the stronger geometric “coupling” 
between pairs of protons. 


V. PERSPECTIVE 


(b) Typical boat-form configuration found only in Ih. 


FIG. 10. The (arbitrarily picked) central proton is labelled 
as “0^^”. The nearest- to fourth-nearest-neighbours are la¬ 
belled “1®^” through “4^^”. From the point of view of the 
central proton, the fourth-nearest-neighbour protons first lead 
to a distinction between chair- and boat-form hexamers. The 
fourth-nearest-neighbour protons are found at distances from 
the central proton ranging from around 3.5 to 5.0 A. The 
green and blue planes are spanned by the relevant 0-0 axes 
connecting the fourth-nearest-neighbour protons. 

populated. This explains why on average the protons in 
Ic move less with respect to their equilibrium positions 
than in Ih, while moving by just as much with respect to 
each other, resulting in the same Ghar as in Ih. 

Going beyond the harmonic approximation, anhar¬ 
monicity reduces the RMS vibrational amplitudes in Ih 
and Ic by around 1% and 2.5%, respectively, localising 
the nuclear wavefunction more in ice Ic than Ih. On the 
level of the collective vibrational modes, the localisation 


Ih and Ic are important in various branches of science. 
Examples include climate modelling and the simulation 
of ice nucleation and formation, where cubic ice plays 
an important role and for which AGanh is an essential 
input parameter. As an example of the relevance in bi¬ 
ological sciences, the benign shape of cubic ice crystals 
is of potential interest for cryopreservation [56| . Here we 
have demonstrated that accounting for anharmonic nu¬ 
clear vibrations is central to understanding and correctly 
predicting the free energy difference between Ih and Ic. 
However, the importance of anharmonic vibrations in 
hydrogen-bonded systems reaches far beyond ice. An 
accurate treatment of anharmonicity is likely to be cru¬ 
cial in correctly describing the energy differences between 
very similar polymorphs of hydrogen-bonded molecular 
crystals which are important in, e.g., pharmaceutical ma¬ 
terials science. 

Calculating anharmonic vibrational energies in solids is 
a challenging computational task, which has only recently 
been successfully achieved using first-principles quantum 
mechanical methods. Anharmonic effects are particularly 
important for light elements, such as the hydrogen atoms 
in H 2 O. Anharmonic vibrations are also expected to be 
important at the surfaces of ice, and when impurities or 
other defects are present. 
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